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INTRODUCTION TO VOLUME II 


This volume consists of a series of appendices whose contents could 
not be included in the body of the report without breaking the continuity. 
The appendices describe topics that were studied during the course of the 
current contract but are only tangentially related to the theme of Volume 
I. Several references were made in Volume 1 to appendices in this volume 
for the details of specific methods. 

The notation introduced in each appendix applies only to that append!^ 
unless otherwise noted. 
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APPENDIX A 
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least squares solution to U.6) will give biased estimates. 
When we use leant squares to minimise che residuals by 

rf- X u * 2 ‘° • 

m K^O 

ve obtain the final expression 


(A.7) 


(A. 8) 


N-l Y-l Y-l 

plo “p Jo ?P+K ‘ " Jo Yn+k 


We will get biased estimates for the a . 

P 


m*0,...,N-i 


(A.9) 



One way to correct this Is to use a method known as Iterative 
generalized least squares. Rewrite (A.7) as 


X a _ m * K"0» 1»««• »Y“1 « 

pfo P P+ K K 


(A. 8) 


Now define new notation for the sake of convenience. First Introduce the shift 
operator q defined so that 


qf ■ f 
q K K+l 


q f m f 

q K K+2 


(A. 9) 


The polynomial operator A(q) Is defined as 


N 


A(q) - j a q® 
m-0 m 


(A.10) 
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Now (A.8) can be written as 


A(q) Y r - W R , let a - a . (A.11) 

In matrix form (A.8) is 


“Y„ Y, Y„1 


~a7 


rw,i 

0 12 


Q 


0 

Y, Y. Y„ 


a. 


W, 

12 3 


1 


1 




■ 


Y 2 Y 3 Y 4 


°2 


W 2 

Y_ Y, Y_ 





L 3 4 5j 


.. _ 




and (A. 11) is 


0 1„ 2. 




r - 

< Y o < Y o ’ Y o 


“o 


w o 

A 


“l 


W 1 




m 


q ° Y 2 qly 2 q2 ^2 


°2 


W 2 

_’° Y 3 «**J 






Now let us operate on W R with the polynomial filter 

B(q) \ m \ (A.12) 

to give n R which is an uncorrelated noise sequence (random variables) so that 
B(q) A(q) Y R - . 

Let A commute with B so that 


A(q) B(q) Y R - n R 


(A.13) 
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and define 


\ - B(q) Y x (A.14) 

so that 

A (q) Y r - n R . (A.15) 

Now equation (A.15) has uncorrelated residuals and hence the solution should 
converge to unbiased estimates. 

The iterative procedure is therefore as follows: 

1. Solve equation (A.11) using the normal least squares procedure. 

This gives an estimate of A(q) which can be thought of as the 

A 

ci's of (A.8). Call that estimate A(q). 

A 

2. Substitute A(q) into (A.11) to produce an estimate of the 

A 

residuals W R . 

A A 

3. Use Wg in (A.12) to make a least squares estimate of B(q). 

A 

4. Calculate Y r in (A.14). 

A 

5. Use the Y R to make a new least squares estimate 

6. Continue this procedure until 

T—* 2 

2, W no longer decreases with the next iteration. 

p-0 K 

The question now is: 

A 

How does one obtain B(q) of step 3 above? 

Assume we have used the least squares process once to find the coefficients a . 

P 



A-4 





Next calculate the y residuals W as 

& 


N 

Z a o ^rvfK " * K-0,1,... ,y-l 

p-o p P+K K 


The results of steps 1 and 2 above gives a y diaenHional vector of the re¬ 
siduals w . 

Calculate the autovariance function of the by 
r(u) « , u-0,1.y-1 


where 


1 Y-u 

c(u) - * I (W 1 - W) <W 1+u - W) , u-0,1.y-1 


- 1 i 

W - * t W . 

Y i-1 1 


Hence the autocovariance function r(u) is defined 

Y (« t - 5; <w t+u - 5) 

r(u) - -- , y-0,1,... ,y-l 

V <w 4 - W ) 2 

i-i 1 


We now have a vector R, y long. 


We want to test this R vector for whiteness of the residuals W^, 
The standard deviation of a single autocovariance function estimate is 


..- iw .i Mr'' r« “11 U ~ 1’i t'lVl f' * * •'* 




The 95% confidence limits for a single autoccvariance sample are 


r(u) + i,96 

ST 

1 96 

Now supposedly if 95% of the samples r(u) are less than ■= — then 

Sy 

we have 95% confidence that the noise is white. 

If the r(u) flunks the whiteness test which will happen the first few times 
we Iterate we must then go back and whiten the residuals WL. 

From Equation (A.12) we want 

B(,) w r * \ • 

L 

Remembering that B(q) - J b q m , 

m«0 



As an example, assume L ■ 2 and K * 3. Then we get 






W 0 w i w 2 

W 1 W 2 W 3 

w 2 W 3 W 4 

Wo W. W, 



Solving Wb - n for the b Q by using least squares gives 


T T 

W A Wb - W n 


b - (W T W) W Y n . 


Therefore the new Y_ previously defined as 




can be written as 


h ■ Vk + Vk+i + V*« ' + V: 


K+L 






APPENDIX B 

EFFECT OF IN^KEASING MODEL ORDER IN PRONY'S METHOD 


A study of the effect of increasing model order in Prony's method 
has been made and the results are presented here along with aome preliminary 
conclusions. In past investigations of Prony's method it has beeu noted 
that increasing the model order above the known order of the waveform 
improves Prony's method's ability to estimate accurately the true poles. 
Although the accuracy of the pole estimates is improved, a side effect 
of this procedure is the problem of distinguishing between true poles 
and those poles that simply fit to the noise and have no relation to the 
Information that is to be extracted from the waveform. 

In this study we relate the inaccuracy or bias of the parameter 
estimates to the amplitude of the residuals of the least-square Prony 
procedure. We assume that the inhomogenous solution (defined in Volume 
I, Section 2) is used to find a parameter vector. The term "residuals" 
is Identical to "equation error". In this appendix, N is the number of 
poles modeled and M is the number of samples. When the residuals are 
zero the least-squares Prony's method either has reduced to curve-fitting 
Prony's method (M-2N) or is processing a noise-free waveform at the proper 
model order. In this case, the pole estimates are quite accurate and are 
free of bias. Conversely the bias in Prony’s method is directly related 
to the magnitude of i residuals. 

Tables B.l through B.9 display the effect of increasing model order 

at three different noise levels. is the standard deviation of the 

N 

noise, is the average standard deviation of the residuals over ten 
Monte Carlo runs. Each table shows the true parameter values in the first 
column, the average parameter values over ten Mone Carlo runs, and the 
variance of each parameter in the third column. For all nine cases the 
time window size is kept constant. Thus, M ■ 100 and AT ■ 0.13 seconds. 


B-l 

















The chree pole pairs shown under the average coition were obtained by 
searching all N poles to find the pairs that were closest to the true 
poles In value. In general, of the extra poles not reported here, N-6, 
have residues two orders of magnitude below those of the true poles. 

From the results of Tables B,1 through B.9, the following observations 
can be made: 

1. Prony’s method seems to guarantee good results If o R < a^. 

2. a n < a., seems to occur when 4N > M. 

That good results begin at o R ■ c?^ Is not surprising. In this case one 
would expect that the residuals are beginning to approximate the noise and 
that the "modeled" waveform approaches the uncorrupted waveform. But the 
observation that c R ■ approximately when 4 n ■ M has no obvious explan¬ 
ation. This observation would obviously not hold true if we were to apply 
Prony’s method to a noise-free waveform. In this case, the residuals 
would drop to zero as soon as the model order reached or exceeded the order 
of the waveform. We might surmise than that thir occurranca at 4N - M 
is attributable to the nature of the uncorrelated noise. 

There is no simple explanation for improved accuracy at higher 
orders. At least two factors seem to be involved: 

1. Increasing the length of the parameter vector (by Increasing 
the order) is an effective treatment for the dense sampling 
problem (which is described in Volume I, Section 3 of this 
report. ) 

2. Making the data matrix more nearly square must reduce the 
magnitude of the residuals. Smaller residuals Implies smaller 
bias in the parameters. 
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Table 3.1. 


Model order study M-100, N-12, o -. 
5 R -.18O0. N 


100 , 


TRUE AVERAGE 



-8.200000E-02 ‘-6. 8 129081-01-7*946325r-03' 

>8*200000E-02 -6. 812908E-01 7.546325E-03 

-1.470000E-.01 -8* 642028E-Q1 9.821274E-03 

-1.470000E-01 ' -6* 64202 8E-Q1-3.821274T*03‘ 

•1*880000E-01 -3* 329218E-01 2.371370E-03 

•1.880000E-01 -3* 329216E -01 2.371370E-03 


“9726 0 0 00 E-01 0* 1 -07- 

•9.26OOQOE-Ol 0* 0* 

2.8740Q0E+0Q 2. 341076E«-00 3.239736E-03 

-2.874000E♦ 00 -2. 3410768♦ 00-STHWCTCTJT 

4.839000E400 4.844311E*00 6.047468E-04 

-4.835000E+00 *4 . 8443UE»00 6.047468E- 04 


' l.OOOCOOEfOO.1. 1093I9E«-00-779'0'9roE=TT2' 

1* QQOOOOE+OO 1*1Q9399E+00 7.509193E-02 

l.QQOOOOE+OO 1*14437 8E + 00 1*381987E-02 

*Ti 0 0 0 00 0 E* 0 0-1.T4437 8EVOO-lT3’mT7E=^^^• 

1* QOQOOOE4 00 1* 189903E+00 6.474729E-03 

1* OOQOOOE+QO 1, 189903E»00 6.474729E-03 


-4*4W111E-“15-47077TOE»2B- 

0* 4* 419111E-15 4* 0 776 92E-28 

0* 1.907967E-01 1.508953E-02 

0*' -1* 5 07967E-01 f75W«E^02 


0* -1* 701996E-01 9*476887E-03 

0. .. 1. 701996E-01 9*4 7 68 87E-03 


REAL PARTS 
OF POLES 


IMAG. PARTS 
OF POLES 


MAGNITUDES 
OF RESIDUES 


RADIAN PHASE 
OF RESIDUES 


B-3 










Table B.2. 


Model order study M«100» N«20, o N «.100 v 
a R ".l24. 


TRUE 

'T.200000E-U2 
•8.200000E-02 
■l*4700Q0E»Ql 
>1*470000£*01 
■1.88OOOOE-Ol 
■1.880000E-01 


AVERAGE 

■9.069350E-02 
•9.069330E-02 
■1* S7178JS?*0JL_ 
•1* f 71783E-31 
-1.951637E-01 
■1* 9916 376-01 


6* 905491E-05 
6.903491E-05 

_1 • 0.BJ3JL 98. E HJ k 

1.083153E-04 
1.278673E-04 
1.278673E-04 


REAL PARTS 
07 POLES 


"*9*2ro00 6 i -01 289997i-01 

-9.260Q00E-01 -9.289997E-01 

2.874000E* 00 2* 883281E+00 

•2.S74800E400 “-2.8832ff 2 «T<f 
4.83900QE+00 4.8402S1E*00 

-4.833000E+00 -4.840211E*00 


3.637123E-05 
3.657123E-03 
1* 0110 44E-0V 

TTiTfnOT=Tf 

1.176039E-04 

1 .1760 39C-04 


IMAG. PARTS 
OF POLES 


l. 00 0 0doE♦'00 ' “' 7 U 0‘27859E♦ 00 2. 0 37443E- 03 

1* OOOOOOE + OO 1* 0 278S9E + 00 2.037443E-03 

„l*m0Q0 E too._l*0.?Afk&71±Jlll_2iliaiA9 fcfli 

1.0 0QQQ0E + 00 1* 0364S7E + 00 2.199189E-03 

l.OOQOOOE+QO 1* 018026E + 00 1.059788E-03 

1.QOQOOOE+QO_ u 018 028 E* DO 1. 0397_88 E-03 


MAGNITUDES 
OF RESIDUES 


O'* 

-1.230493E-02 

9.62Q293E-04 

0* 

1.230493E-02 

9.820293E-04 

0* 

-1.928778E-02 

2.042011E-03 

0* 

l.^SttOE^OE 

2.042011E-03 

0* 

-E.326322E-03 

2. Q243 03E-Q3 

0* 

5.328322E-03 

2.024303E-03 


RADIAN PHASE 
OF RESIDUES 


B-4 
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Table B.3. Model order study M-100, N-36, o„-.100, 
3„-.0678. 


TRUE 

8.200000E-02 
6 * 2000002*02 
1.470000E-01 
1.470000E-01 
1.8800002-01 
1.8800002-01 


AVERAGE Q 

-6*. 169175*2-02 3. 0566572-05 

•6.1691712-02 3.0566572-03 

.. -1.>5353.52-01_VO7.227JE.-0> 

-1.4835382-01 1.1722732-04 

-1.0425902-01 1.8613972-04 

._=0fJ»mUbll_UM 1W1-0 4 


REAL PARTS 
OF POLES 


9. 260 0OOf-Oi" 
9.260000E-01 
2.6740002+00 
2.8740002+00 
4.6350002+00 
4.6350002+00 


•9.2664802-01 2.3949362-03 

2.8748362+00 5.3018692-05 

• 2.8 n a m ♦ on -573*1)18OT2- 0 i 

4.8361282+00 4.4500362-05 

•4.83612SC+0Q 4.4500362-03 


IMAG. PARTS 
07 POLES 


1 . 0000002*00 
1 . 0000002*00 
1 . 0000002+09 
1 . 0000002+00 
1 . 0000002+00 
1 . 0000002+08 


9.9310762-01 1.3250782-03 
9.9310782-01 1.3250782-03 
.1.0137202 +00 4.1 19 6332-0 3 
1.0137202+00 4.1196332-03 
9.6910632-01 6.2767592-03 
9. 6910632-01_6.2767592-03 


MAGNITUDES 
07 RESIDUES 


•3. *5631472-03 
3.5631472-03 
1. 445358 2-03 
• 1.445*582=11 
•1.4026732-02 
1. 4026732-02 


8.5535592-04 

8.5535592-04 

IjilffttibJUL 

177902252-03 

7.1345402-04 

7.1345402-04 


RADIAN PHASE 
07 RESIDUES 













Tfebla B.4. Modal ordar itudy M-100, N-8, c^-.00100, 
o R -.007098. 


TRUE 

a NBA. flUA -A ~T M AM » IBM 

AVERAGE 

a 

tnim-TT* 
-8.2000008-02 
-1.4700006-01 

-5774682)8-01 
-4, 7460212-01 
-3.23707)8-01 

^.f7l5TO?33 
6. 6719 06E-03 
1.1543838-03 

•1. 47Q0fltJE-()i 

-l.seooooE-oi 

-1.880000E-01 

-3.217'07nr“ffT" 

-1.9949998-01 

-1.9949191-01 

“Tm'4irnnr 

2.614657E-P5 
2.614697E-0) 

• 

9.2600008-01 
-9.260000E-01 
2*874000T4 0-0. 

8.99030)8-01 

-6.99030)8-01 

2.8743948+00 

8.704512E-04 
6.704912E-04 
2.6016898-04 

-i.iHoooeioo 
4* 8390008400 
-4.8350008400 

-2.874394I+00 
4.6432412400 
•4.8432418400 

2. 6016892-04 

3. f698662-0) 
3.669866E-0) 


1« 0uuOOu c □ 
1.000000E+QO 
1.0000002400 

177371898400 

1.7371998400 

1.2238202400 

9.114346E-0T 
9.1143668-03 
1.4617208-03 

l* ooooooEioil 
It OOOOOOE + OQ 
1.OQOOOOE+QO 

raimrur 

9.0263708-01 

9. 02837 0e-01 

1.4)l720K-03 
4.5900 50E-04 
4.950050E-04 


REAL PARTS 
OF POLES 


IMAG. PARTS 
OF POLES 


MAGNITUDES 
OF RESIDUES 


"0.- 

2. 0 228 031-01"“ 

“77ZT5594T-03 

0. 

-2.0228038-01 

7.2359948-03 

0. 

2.29633)2-01 

3.769424E-03 

0. 

-2.2983j)X-01 

2.7694248-03 

0. 

1.2808948-01 

2.4557898-05 

0. 

-1.2606948-01 

2.49576IE-0S 


RADIAN PHASE 
OF RESIDUES 



Table B.5. Model order study M-1G0, N-12, a 


TRUE 

3 r -. 00300. 

AVERAGE 

A 

-8.200000E-02 

-6.3347606-02 

3.S19636E-Q5 

•6.200000E-02 

•6.3347S0E-Q 2 

3.819636E-06 

-.1 *.47.0 0 001-01_-11.482 0 4i E-.0i_. 

_3*.3.0 75561-06. 

-1.470000E-01 

-1. 482045 E-01 

3.30 75 56E-QS 

•1.66 OOOOE-Ol 

-t. 679567 E-01 

1.655490E-07 

| 

I 

■ 


,.JkjJUU&9UL6±0JL 


9.260000E-01 9.2632276-01 

-9.260000E-01 -9.2632276-01 

if* 835000E+00 4. 6390411+00 

-4* 6390Q0E+Q0 -4.S3904II+00 


3.0542061-06 
3.0 542 06E-06 
2.31 3 3 36E-0 6 

iltimaHi 

6.7570916-07 
6* 7570516-07 


1.OOOOOOE+OO 1.00544L6+00 9.2241066-05 

1.00 Q0QQE + 0Q 1.0054416+00 9.2241066-05 

^l.aMflJaE tPQ ., 4.7-160 98W1 

1.OOOOOOE+OO 1.0038536+00 4.7160966-05 

1.OOOOOOE+OO 9.994317E-01 4.064030E-06 

l.OQOO Qflg+QO .9 . 994 3 171-01 4. 06 40 306-0 Bl 


0* 

-1.1104116-03 

6.015323E-05 

0. 

1.11041}6-03 

8.0153236-05 

0. 

6. 180 2731-04 

4.066317E-0S 

0. 

-8.1802706-04 

4.066317E-05 

0. 

1.6991516-03 

9.113711E-06 

_1*_ 

-1. lilUIMl 

... 9.1137111-06, 
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-. 00100 , 


REAL PARTS 
OF POLES 


IMAG. PARTS 
OF POLES 


MAGNITUDES 
OF RESIDUES 


RADIAN PHASE 
OF RESIDUES 





Table B.6. Model order study M-1G0, N-20, a N -.001, 
cr R -. 001246. 


TRUE 

•8.200000E-02 
*8*20 00006-02 


AVERAGE 

'8. 199893E-Q2 
•8. 199895E-02 


6.9594356-09 

6.939433E-09 


wuww * 

-1.470000E-01 
-1.88OOO0E-O1 
■li98 0QQ0E-Q1 

V. Zl ***** y * 

-1.47015SE-01 

-1.8303S8E-01 

-1.880399E-01 

1.0346476-05 

1.1212236-08 


9.2600Q0E-01 
-9.260000E-01 
2.8740QQ£*00 

9.2600796-01 

-9.2600798-01 

2* B7402IE*00 

3.504264E-03 

3.304264E-03 

9.5360526-09 

-2.8740006*00 
4* 835000E*00 
-4.-833Q_0P_L»M_ 

-2. 874025H-00 

4. 8349946*00 
-4.83499*6*00 

9.536032E-09 
1.374833E-08 
1« 3745336-09 


1.0000006*03 
1.0000006*00 
1* 0000006*00 

9.9995096-01 

9* 9995396-01 
1.0001106*00 

2*1308856-07 
2.1308536-07 
_ 2.J.91616E-Q7 

i.OOQOOOE*OQ 
1* QOOOOOE>QO 
1. 0000006*00 

1* 0001106*00 

1.0000666*00 

1. 00006Se*00 

2.0916146-07 
9* 0947372*08 
9.0947376-08 


0. 

0. 

-4* 4011026-05 

4* 4019022*09 
-8.9732606-05 

9.Q35431E-T8 
9.0354316-08 
2. Q 045726-07 

0. 

0. 

JU _ 

6.9732906-05 
-2.9039196-03 
2.9Q39UE-05- 

2. 0049726-07 
2.2595966-07 
-.Jiili ZSS&zAL 


REAL FARTS 
OF POLES 


XMAG. PARTS 
OF POLES 


MAGNITUDES 
OF RESIDUES 


RADIAN PHASE 
OF RESIDUES 


B-8 



Table B.7. Model order study M-100, N-8, O..-.0100, 
o r -.0394. 


TRUE AVERAGE a 2 

8/20 0000£-02 -4, 81594 IE-01 4.313S74E-Q3 

8« 20 Q000E-Q2 -4.815941E-01 4.313874E-03 

1.470000E-01 -4«931019 E-Ql 1.633 2 14E-03 

1.470000E-01 -4.931013E-Q1 1.6332UE-03 

1.8800005-01 -2* 503945E-Q1 4.518706E-05 



REAL PARTS 
OF POLES 


9.2600005-01 
9.260000E-01 
2.87400Qg»0Q 
2.874000E4-00 
4*83S000E+ 00 


■2.44487SE+0O 

4.820967E+00 


6.292900E-04 
1.699293E-04 


IMAG. PARTS 
OF POLES 


l.OOOOQOS^OO 8.254609E-01 2.930682E-Q2 
1* Q00080E+ 00 8*254609E-01 2.930682E-02 
1.00 OOOQE»OQ 1*154Q58E+Q0 3.968307E-03 
1 . 000000**00 l.tiVoV 8 E 400 3.S58307E-Q3 
l.OOOOOOEfOO li112501E^00 4.643365E-Q5 
1»0000005*00 1.112301E+00 4.643363E-05 



MAGNITUDES 
OF RESIDUES 


2* 991933E-15 
1» 623 083 E-01 
•1.6230S3E-01 
•1. 075143 E-01 


6* 5 963*»lE-28 
2. 7298 86E-03 
2.729888E-03 
9. 082337E-04 


RADIAN PHASE 
OF RESIDUES 


'l 


I 


I 


ia 
















T^ble B.8. Model order study M-100, N-12, g «.0100, 
3 0295. 


TRUE AVERAGE g 2 

8.2 0 0 0 0 0 E-02 -2.35316iE-fll 2a <3592 68E-03 

8. 2000005-02 -2.35315 IE-01 2.859268E-03 
1.470POPE-01 -2. 4074SSE-Q1 5.513517E-QV 

1. 470000E-01 -2.407455E-Q1 5.518517E-0'. 

l.anOOOOE-Ql -1.9345&1S-01 1.094043E-05 

1.SHO000E-Q1 -1.934551E-01 1.094Q43E-Q5 


REAL PARTS 
OF POLES 


9.26 0 00 0S-01 
2.874000E+Q0 

9.4d447J£-di 
-9.408873E-01 
2. 899695 E+ 00 

JTT5TT53^ir 

3.593158E-04 

1.687968E-04 

^2737^ 000 ETO fl~ 
4.835000E+00 
•4,8390005+00 

rnnr 

4,842363 E+00 
-4,842353E+00 

r." , 6'87g , 6 , 4i-u; 
1,122618E-Q4 
1.122S18E-Q4 


1. Q0OQ0QE+QQ 
1.000000E+00 
1* OOOOOOE+OQ 

1.4i60SQE+00 
1* 426050E+0 0 
1.211075E+00 

1.70 80 38E-02 
v 1.70803BE-02 
v* 2,967772E-03 

1, OOOOOOE+OQ 
1.OOOOOOE+QO 
1.OOOOOOE+OQ 

1.211075E+00 
9.755319E-01 
9.795319E-01 

2.967772E-03 
4,1906 53E-04 
4.190653E-04 


s. nmye-ai — Tttzrunn 

2. 712845E-03 7.626122E-03 
7* 308093E-Q2 3.766332E-03 

77TJnTT5T^ff?- 3; 7 6 T 332E-d3 

1.0259S3E-Q1 8.128753E-0V 

1. 025580E-01 8*1287 53E-04 



IMAG. PARTS 
OF POLES 


MAGNITUDE 
OF RESIDUES 


RADIAN PHASE 
OF RESIDUES 




1 


i 

I 
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APPENDIX C 

THE ADAPTIVE METHOD FOR RESONANC E ESTIMATIO N 

The adaptive method reeulta from the generalized scheme of Volume I* 
Section 2, under the following assumptions: 

1. a “ 1 and g *0. 

o o 

2. Pj ■ —-— , i ■ l,...n and P ■ 1. 

i z-z^ * o 

3. The model input is a unit sample at k - 0 (discrete impulse). 


i 

} 


The unique feature of the method is that the filter poles, z^, may be 
adjusted to any value in the Z-plane, An adaptive technique fnr adjusting 
the filters consists of first initializing the filters to arbitrary values 
in the Z-plane and repeating the following steps: 

1. Find an estimate of the process transfer function using 
the current filters in the model. 

2. Set each filter pole to one pole of the estimated 
transfer function. 

This procedure is repeated until a 1 approach zero. The poles of the 
process can be estimated during the course of iteration by 

A *i 

z i 1 + a t 


removing the need for finding the roots of a polynomial. The filter poles 

are updated to on each iteration. When the a^ approach zero the pole 

updating ceases and the method converges. At each iteration, the a. are 

1 

found using any of the techniques for finding a parameter vector described 
in Volume I, Section 2. Perhaps the simplest method is to choose the 
parameter vector as the weakest eigenvector nf ft*Q. At convert jnce the 
S-plane poles, can be obtained from the fi’.tsr poles by 


s 


i 


in z^ 
T 


and the A^ • 3^ 
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The adaptive filtering method seems to be a very useful technique 
because it uses what seems to be optimal filters, band pass filters, for 
estimating pole locations. In addition it presents a measure of the error 
and iteratively adapts the filters until the error is at the level where 
it is desired. 

As an example of the use of the adaptive filtering method consider 
the data shown in Figure C.l. These data were numerically generated by 
using the time domain computer code TWTD [C.l], The structure modeled 
by TwTD was a thin cylindrical scatterer. The signal was contaminated 
with noise to give a 15 dB signal-to-noise ratio. Figure C.2 shows the 
resulting poles from five Monte Carlo trials. 

The following statements can be made about the adaptive method after 
studying it. 

The adaptive method is a new method which, in many cases, pro' ides 
excellent pole estimates under difficult conditions. The method is 
unique in that a solution to a polynomial is not required to find estimates 
of the process p.les. the method, in effect, "swallows" the polynomial 
solver in its own iterative pole-searching scheme. 

Attempts to analyze waveforms consisting of highly damped exponential 
components, such as the transient responses of a sphere, have not been 
successful. The adaptive method does not converge for waveforms which 
display double pole characteristics, that is, waveforms with components 
of the form t exp(st). Slight modifications to the adaptive method might 
allow the analysis of such waveforms. 


C-2 





































References 

[C.l] M.L. Van Blaricum, "A Numerical Technique for Che Time-DependenC 
Solution of Thin-Wire Structure with Multiple Junctions", 
Electromagnetics Laboratory Repo rt, 73-15, University of Illinois 
Urbana, Illinois, December 1973. 







APPENDIX D 


THE PENCIL-OF-FUNCTIONS METHOD 

Tha pencil-of-functions msthod rssults from ths generalized modsl 
describsd in Volums I, Ssction 2 with ths following assumptions: 

1. Oq ■ 1 and 8 q ■ 0. 

2. ■ F^Cs) ■ (l/s) 1 (cascadad continuous intagrators). 

This mathod is not vary assy to implamant on a digital computar sinca tha 
continuous-tirna intagrators cannot ba implamented exactly by any algorithm. 
Whan approximate integrators ara cascaded, as they are in tha pancil-of- 
functions mathod, largo errors can ba quickly accomulated and tha intended 
result after a number of integrations destroyed. 

This difficulty can ba resolved by cascading discrete intagrators to 
obtain filters with pulsa transfer functions given by: 

r t * 'i (,> ■ (A) ■ (if 

which can ba implemented on a digital computer with no error by using 
difference equations. The variable Z is defined as 


2 " 1 * 


and z is the z-transform variable. 

Whan discrete integrators are used the discrete pencil-of-functions 
mathod results. The poles of the pulsa transfer function of the process 
or waveform may be estimated as 


s i “ l-V 















kL 

where is ths i zsro of 

* n + a. Z n-1 + ... + a - Q (D.l) 

l n 

In ths original pencil-of-functions mtthod ths wars found by using 
Jsln's msthod for constructing ths parameter vactor which is dascrlbad in 
Voluma I, Saction 2. With Jain's mathod 



aU ak 

whara A.., in this csss, is tha olsmant at tha i row and j column of 
»dj fl*n. Othar methods can ba usad to astimata tha paramatar vactor. 

Tha S-plane astimataa of tha poles, s^, ara ralatad to tha zaroa of 
(D.l) by 

i 

Ona difficulty with tha pencil-of-functions mathod is ralatad to 
tha attenuation of tha higher frequency modes of tha process output by 
tha rapaatad integrations applied to the output waveform. It can ba 
verified that an integrator is simply a first order filter whose Laplace 
transfer function has a pole at tha origin in tha S-plana. Such a filter 
tends to suppress tha higher frequencies present at its input. Tha 
higher frequency suppression phenomenon is illustrated in Figure D.l. 
Normally, when an exponential function is integrated repeatedly, 
components of power of time exist in the higher integrals as wall as the 
original exponential function components. In Figure D.4 the components of 
powers of time have bean subtracted from tha integrated waveforms in 
order to make the attanuation of the higher modes more evident. Tha 
first waveform is a hypothetical waveform provided for analysis. Tha 
waveforms that follow are the integrals of Increasing order of tha first 
waveform and display the increasing dominance of the fundamental mode or 
mode of lowest frequency. Further integrations yield nearly identical 
waveforms. The integrated waveforms tend to become linearly dependent 
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vavaforma. Tha intagratad wavaforma tand to bacoma linaarly dapandant at 
highar nodal ordara. Tha matrix than tanda to alngularity and tha 
mathod bacomaa unatabla. Tha auppraaalon phanomanon occura in both tha 
dlacrata and continuous mathoda. In fact, avan for vary modaat modal 
ordara, tha mathod can bacoma numarlcally ill-conditionad to a dagraa that 
apacial cara muat ba takan to aaaura accurata invarslon of n*n. 









APPENDIX E 


METHOD OF REDUNDANT AVERAGING 


Tha redundant avaraglni schema la a preprocessing schema Chat attampta 
Co oomblna redundant data (that la, more data than ara necasaary to determine 
tha paranatara) In a way that avolda tha blaa lntroducad by ualng a laaat- 
aquaraa achama to comblna radundant data. Tha mathod attampta to tranaform 
a raw wavaform of aora than 2N aamplaa whara N la tha modal ordar Into a 
praprocaaaad wavaform of exactly 2N aamplaa by avaraglng within tha wavaform. 
Tha avaraglng can ba dona ao that tha axpactatlona of tha polaa of tha pra¬ 
procaaaad wavaform ara aqulvalant to tha axpactatlona of tha polaa of tha 
raw wavaform provldad tha addltlva nolaa on tha raw wavaform la aaro maan 
and uncorralatad batwaan sucasalva aamplaa. Tha praprocaaaad wavaform may 
than ba procaaaad with cuxva-flttlng Prony'a mathod to avoid tha blaa of tha 
laaat-aquaraa procadura. 

Tha moat ganaral daacrlptlon of tha radundant avaraglng achama can 
ba atatad almply aa 

v 1 

' * jlo Vu<, ’ k " 0 ’ 1 . 21,-1 

aU 

whara N la tha modal ordar, X. danotaa tha k aampla of tha praprocaaaad 

wavaform, and y^ danotaa tha k aampla of tha raw wavaform. In ordar to 

limit thla daacrlptlon to tha aaaantlal faaturaa of tha radundant avaraglng 

achama, N., tha numbar of avaragaa, an'' N , tha dacimatlon apoch, ara not 
A ■ 

axpllcltly dafinad hara. Thaaa parameters ara chooaan aa daalrad but uaually 
in a way that would produca a daalrabla praprocaaalng moda in soma sanaa. 

For inatanca, tha value for tha dacimatlon epoch might ba choaon on tha 
basis of tha maximum frequency, <u , praaant in tha raw wavaform (if thla 
information is available) and tha numbar of avaragaa might be choaan ao 
that every aampla in tha raw wavaform la used once tha value of tha 
dacimatlon apoch la sat. Hanes tha sampling interval, At, for tha 
praprocaaaad wavaform could ba obtained aa 

















irfi’iir* 'iiiuaa^t 


At - 

where ■ 

The number 


N At 

• raw 


Integer Part 



of avaragaa, N^, can ba computed by 


N A - M - N- (2N-1) 


whara H la tha number of aamplea In the raw waveform and N la the order 
daalred. It ahould ba noted that tha above method for determining and 
Nj la not tha only poaalbla technique. 


A 






Tha redundant-averaging procedure la affectively two operational 

1. Apply a low-paaa moving-average filter, A(z), to tha raw waveform. 

2. Decimate tha filtered waveform with decimation epoch N g . 


Tha tranafar function of tha moving-average filter iat 

N.-l N.-2 K i 

A(a) - a * + a A + ... + a+1 • 


The numarator of A(r) can be factored into 


II (z-z.), 
i-1 1 


The zeroa, z^ of the numerator of A(z) have unit magnlcudaa and argumanta, 


erg z t 



for i ■ 1, N^. It la clear that tha firat factor cancala with tha 

denominator giving 


I 


<1 


< 

3 

1 

i 


Na 

A(z) ■ H (z-z.), 
i-2 1 
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ZEROS OF A(Z) 


L ' 

Figure E-l. Z-plane plot of 

j 

preprocessing f 
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l' ** 


r 



M 


























Table E.l True poles and results of five Monte Carlo runs for the 
redundant averaging example shown in Figure E.l. 


REAL PARTS OF POLES 


TRUE 

POLES 


MONTE CARLO RUNS 


2 ' 

3 

-.127 

-.122 

-.157 

-.239 

-.103 

-.227 


AVER. 


-.107 

-.187 

-.271 


Z DEV. 


IMAGINARY PARTS OF POLES 


MONTE CARLO RUNS 


TRUE 

POLES 


.926 

2.874 

4.835 


2 


.904 

2.913 

4.899 


• .942 
2.889 
5.065 


5 


.936 

2.882 

4.803 


AVER. 


.925 

2.902 

4.899 


X DEV. 








































A(z) chen has zeros which are evenly spaced around the unit circle but 
the zero at z»l is canceled. Since there is never a zero at z»l, the 
moving-average filter can be viewed as a type of low-pass filter, A 
z-plane plot of the zeros is shown in Figure E.l. 

As an example of the application of this method consider the noise- 
free waveform shown in Figure E.2a. This waveform consists of 400 data 
points representing a sixth order exponential function generated with the 
poles listed in Table E.l. This data was then corrupted by adding white 
noise with a standard deviation, o, of 0.5. Figure E.2b shows the noise 
contaminated signal. This signal has a signal to noise ratio of 15.6 dB 
where the signal to noise ratio is defined as 

S/N - 20 log 

where R is the peak amplitude of the transient signal, Five Monte 
Carlo trials were run on this data using the redundant averaging mathod. 

The model order was selected to be 8 (two more than the known order) and 
Ng and were chosen to be 16 and 160 respectively. Figure 3.2b shows 
the reconstructed waveform obtained from one of the Monte Carlo trials. 

Note how good the fit is considering the signal to noise ratio was 15.6 dB. 
Table E.l lists the results of the five trials and shows the pole averages, 
standard deviations and per cent deviations. It should also be noted that 
the signal to noise ratio compared to each pole residue Is only 6 dB. 

Hence from poles with a 6 dB information content we were able to recover 
them very accurately. 

The redundant averaging scheme has two difficulties or limitations 
that can cause the method to be lass effective: 

1. If the sampling rate in the preprocessed waveform is 
sufficiently low, the higher frequency poles can be 
folded, perhaps several times, About the Nyquist 
frequency. 
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2. The method can null modes in the preprocessed waveform 
that exists in the raw waveform. 

The frequency folding phenomenon introduces an ambiguity in that it 

is not known how many times a particular pole has been folded about the 

Nyquist frequency. Hence, N possible poles are introduced for each 

s 

extracted pole by the redundant averaging scheme where 

N » 

3 


(AT) preprocessed 
(AT) raw 


Of course, if the highest frequency mode of the waveform is known to be 
lower in frequency than the preprocessed waveform's Nyquist frequency, the 
ambiguity is resolved but, in general, this will not be the case. 

The mode nulling can occur if the averaging parameters are such that 
the zeros of the resulting low pass fi-ter cancel poles of the signal 
Itself. That is, if we define the preprocessed waveform P(z) as 

M 

P(z) ■ A(z) D(z) 

resulting from the averaging process A(z) operating on the signal D(z), 
then if A(z) and D(z) each have terms in common they can tend to cancel 
each other. 
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COLUMN PRONY'S METHOD 

la this method the z-plane estimates of the poles are the roots of 
the polynomial 

„ , N. , , N. 2 . > N. N . 

o'q + ct^Cz ) + a 2 (z ) + ... + a^U ) - 0 

The are determined from the system, 



where ct^ is assumed to be one and {y^, y^, y^ 2 + jj_ 2 _} denotes the 

sequence for the waveform containing N^+N samples. It should also be 
noted that the. polynomial has N^ roots only N of which are related to 
estimates of the true constituent poles of the waveform. 

2 

Column Prony's method offers n means of combining N +N data points 
compared to the 2N data points of the standard Prony's method. The 
column Prony's method can either be used in the least-squares version or 
the curve-fitting version. If the curve-fitting version of column Prony's 

method is used, the resulting parameter estimates are unbiased. 

2 

Unfortunately the method yields N pole estimates only N of which are the 
true poles. It appears then that column Prony's method leads to the same 
problem that increasing the model order led to in the standard Prony's 
method: an ambiguity in the identity of the true poles. Moreover, this 
method has an additional problem: the matrix can become nearly singular 


P-1 


» 








even if the model order is lower than or equal to true order. This 
phenomenon is similar to the singularity of the standard Prony's method 
when the model order equals the waveform order and the highest frequency 
is equal or nearly equal to the Nyquist frequency. 
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EVAN'S AND FISCHL'S METHOD 


I 


» 


» 


► 


» 


In Evan's and Fischl'a method [G.l] they define a "true error" 
sequence {e^, k*0, .M-l} as the error between the given waveform and 
the "fitted" waveform. They then proceed to define the "equation error" 

sequence {d^, k*0.M-l}. They further proceed to define a relation 

between the equation error and the true error: 

e ■ Wd 


where 


o*i Vi ] 
i 3 td o d i •“ d M-N ] 


and W ■ A[A T A] ^ 


where 


r* a. 


“N a N-l 


UO 


... 0 “1 


• • e • 


0 

a 0 


■ • • I 


a N-l 
°N - 


M is the number of samples, N is the number of poles, a^*l, and the a^; 
are the coefficients of Prony's difference equation defined in Section 2, 
Volume I. 










By ueing this relation between the two errors they describe two iterative 
procedures aimed at minimizing tha "true" error criterion: 






Ic can be shown that 



W(a) - 


3a 


N-k 


A(a) 


- W(a) 


3a 


A(a) 


N-k 


A(a) 


+ A(a) 


3a 


N-k 


Av .) 




A(a) 4 A(a) 


and [(3/3a N-lt ) A(a)] is simply th* matrix A(a) with onaa raplacing tha a N _^ 
with all othar alamants saro. 

Tha following observations can ba mada about tha mathodt 


1. It produces optimal pole estimates in tha sense that it 
minimizes "true error". This means that tha mean-square error between tha 
given wavaform and tha fitted waveform is minimized. These optimal pole 
estimates are obtained only in the second iteration phase. 

T 

2. .ue matrix A A must be inverted on each iteration of both tha 

T 

first and second procedures. When the wavr. . u has over 100 samples, A A 
becomes very large, and hence, expansive tu invert. Consequently, the 
method is prohibitively expensive when the wavaform has over 100 samples. 
This method has yet to be evaluated by tests on noisy data. Nothing is 
presently known about its convergence characteristics of its tolerance of 
noise. 


The following example illustrates the optimal estimates. 
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Exaaola 


N • 1, M - 3 

2" [2 1 2] (vavaforo eo ba approxiaaead) 

s - lju- t"jl 

Initial paramatarai , a ■-!, • Thasa ara cha optimal paramatara, 

tharafora Cha mathod should convarga inanadlaCaly. 



'2/3 1/3' 
.1/3 2/3. 







rrx 

0 


-2/3 -1/3 - 

/ _ 

-1 0" 





/ 

10 0 



0 

l 


1/3 -1/3 


0 10 


1 -1 

IL° 

0_ 


1/3 2/3 

' 

0 1 


» 


+ 


1 1 o' 
0 -1 1 . 


I 


1 0 
0 1 
0 0 


2/3 1/3' 
1/3 2/3 
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Therefore the updated valua of la equivalent to ita old value, hence the 
method hee converged at the optimal parameters for the second Iteration 

phase. 
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APPENDIX H 

CONSTRAINED PRONY'a METHOD 
H.l Uni of a Constrained Prony Method 

Many timaa, In tha application of tha Prony algorithm aoma of tha 
polaa ara known a priori. It would ba vary uaaful if tha Prony algorithm 
could ba conatralnad in auch a mannar that tha knowladga of tha known 
polaa la uaad in axtracting tha unknown polaa. Tha known polaa could ba 
polaa of tha driving function which ara known from knowladga of tha Laplaca 
tranaform of tha driving function, ayatam polaa which ara known from pravloua 
Prony analyala or othar tachniquaa, or polaa introduced to modal tha nolaa. 

It la wall known that for cartaln data aata Prony'a algorithm haa 
aotua difficulty In axtracting tha trua polaa that ara containad in tha 
data. Thla problem ia ganarally ralatad to tha nolaa in tha data but will 
not ba dlacuaaad hara. Any mathod by which tha accuracy of tha trua polaa 
can ba incraaaad will ba vary uaaful. 

Logic would tall ona that making uaa of known information ahould 
incraaaa tha accuracy of a calculation. Hanca, if uaa la mada of known 
polaa ln tha data it would aaam raaaonabla to aaauma that aoma of tha 
lnatabillty ln Prony'a mathod would ba allavlatad. Tha proof of tha 
validity of thla atatamant raata on tha actual lmplamantatlon of tha 
conatralnad Prony algorithm and tha raaulta comparad for tha aama data 
analyzad by tha unconatrainad mathod. 

In addition to aiding ln lncraaaing tha accuracy of tha trua polaa 
it might ba poaaibla to introduca random polaa which will modal tha nolaa. 

In that mannar, tha polaa that ara knoim a priori ara not tha polaa which 
wa aaak. Tha difficulty with thla, if it ahould work, ia that wa naad to 
know tha rank of tha ayatam ao that wa can introduca tha propar numbar of 
nolaa polaa. 












Another asset of having a contrained Prony's method Is that it can be 
used to aid in deconvolution. The following discussion follows directly 
from work performed in Reference [A.2]. 


Assume that a system is excited by a driving function which can be 
represented in the form 


M 

G(t) - 2 
j-1 


8 j 



u(t). 


(H.l) 


It io presumed here that the driving function poles s^ and the associated 
residues g^ are known. These can be determined analytically if the ana¬ 
lytical form of the driver is known or can be determined from a Prony's 
method fit to measurement data of the driving function. 


Now assume that the response of the system to the driving function 
G(t) is 

M+N st 

R(t) - T r e 1 u(t) (H.2) 

i-1 1 

and that the impulse response of the system is 
N s t 

H(e) - I h e k . (H.3) 

k-1 * 

Hence, there are N system poles and residues and M poles in the driving 
function. Of course, the response function R(t) can be written as the con¬ 
volution of the driving function G(t) with the impulse response H(t). That 
li» 

R(t) - H(t)*G(t) , (H.4) 


where * denotes convolution. 





What is usually desired is to obtain tha N system poles and residues 

of (H.3). This must be done by deconvolution of the driving function 

(H.1) from (H.2). If the r^ and s^ are known, the MfN values of Sj, contain 

the M values s.. The g, and the s. are known, then the N values of h and 
J J j k 


3 k can be determined analytically. 


It can be shown that the response function can be written in terms of 
the driving function and the impulse response as 


M N \ ■ t N M g. s 

I ( «i I rrr>• J - I c-t I 7-£->* 

J-l 3 k-i s k k-1 k >1 a J *k 


(H .5) 


Hence, the response function residues can be defined as 

N h k 

r, - 8, 1 • for 1 " 1 • M 

1 J k-1 8 j a k 


(H.6a) 


r, - -h V ——L- , for i - M+l , N+M 

1 * j-l V k 


(H.6b) 


From (H.6b^ the residues of the impulse response are written as 


h k-|3: 

j-i 3 j 


k - 1 , N 


i - M+k 


(H. 7) 


Thus the res' *tes or amplitudes of the impulse response can be obtained from 
the known va. js r^, g^, s^, s^. Of course, since the M values of are 
known and the M+N values of s^ are presumed to contain the M values of , 
then the N values of s^ can be determined by inspection. 









Prony's method, however, will not necessarily give the true values 
of the driving function poles in the extracted poles of the response 
function. Hence a constrained Prony’s method is required so that analyti¬ 
cal deconvolution of (H.7) cau be obtained. 

Work performed so far Implements Method 1 which is described below. 

It was found that this method works perfectly as long as no noise is 
present in the signal. As soon as noise is introduced in the signal and 
a least-squares method is used the matrix (H.17a) is corrupted with noise 
which through the least-squares process also corrupts the constrained 
parameter. Experiments have confirmed this observation. This suggests 
that the constraining Method 1 may work well in noisy data if the curve¬ 
fitting Prony's method is used. This has not been tested as yet. 

Method 2, forcing the polynomial root solver to find certain roots, 
has also not been tested. This is because if the coefficients are all 
corrupted by noise through the least-squares procedure then subtracting 
any given poles out of the polynomial will force the remaining poles to 
carry the burden of all the noise. 

H.2 Method 1 


In the implementation of Prony's method, an Nth order polynomial is 
solved for its roots. The order of the polynomial, N, is the number of 
poles being sought in the transient data. If the coefficients of the 
polynomial are denoted as a then the polynomial can be expressed as per 
Reference H.l as 

a Q ♦ ajZ 1 + a 2 Z 2 + ... + a N Z N - 0 (H.8) 
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where is usually set equal Co unity. The N roots, Z 4 , are defined as 



with s^ being the poles sought, and At being the time step size used in the 
analysis. 


If the value of one or more of the poles is known - that is, we know 
some of the - then the can be substituted into (H.B). For example, 
if s^ or Z^ is known, then (H.B) can be written as 

Og + <*2^i * * * "** a N^i " ® • (H10) 


The N+l polynomial coefficients a are solved for in Prony's method by 
solution of the difference equation 


N-l 


p-0 


P P+K 


-I 


N+K’ 


k-0,l,...,Y“l 


y-m-n 


(HJ1) 


The and I^ (R are the samples of the transient signal being analyzed and 
M is the total number of samples being used. The value of M must be at least 
equal to 2N to give N sets of equations in the N unknowns a p . However, if 
the value of a pole is known, then one of the N aquations can be equation 
H.10 and N-l equations of the form of H.ll can be written. 


If L poles are known a priori, the L equations of the form of (H.10) 
can be written as 


N-l „ „ 

V a z p «-z, £ — 1 , 2 , ,.,,L 

p l l 

p-0 v 


(H.12) 
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and y»M-N-L equations can be written in the form of (H.ll) 

N-l 

p-o ° p Ip+K " " Ik+n * k " 0,;L . Y_1 ’ Y " M ” N “ L • (h.13) 

Hence there are still y-M-N total equations to solve for the N values of a 
however the system is constrained by the knowledge of the location of L poles. 
As is usually done in Prony's method, if M-2N then the set of equations is 
inverted and solved. If M>2N then a pseudo-inverse procedure is used. 


Using the matrix notation of Reference H.l, if M-2N and L poles are known, 
then we solve the equation 


AB - C 


(H. 14) 


where A is a square matrix defined as 

1 Z 1 

1 Z- 


N-l-L 


N-L 


• • » Z 

• • • Z 


N-l 

1 

N-l 

2 


.N-l 


... I 


N-l 


... I 


N 


^-l+l *•* ^2N-L-2 


(H. 15a) 
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and B and C are vectors defined as 


°N-1 


(H.15b) 


(H.15b) 


L 2N-1-L 


If M>2N and L poles are known, then the solution takes on a pseudo-inverse 
or least-squares form as 


T T 

A AB ■ A C 


(H. 16) 


Where A la the transpose of the matrix A and A la now a rectangular matrix 
of the form 

i z : z* ... zj’ 1 

l z z 2 2 n- x 

i i>2 u 2 ••• Z 2 


M-N-l-L 




... I\ 


... I. 


... I 


M-L-2 


(H.17a) 








(H.17b) 


Consider the simple example of a two pole system 

st at 

f(t) - 2e 1 + 3a L 


where s^ ■ 0 and s 2 • -A. Assume that s^ Is known, thus giving Z ^ 
Also, let f(t) be sampled at At - 0.2 seconds, giving a data set as 


2.6057 


2.2722 


2.1223 


Using the constrained Prony's method in the square system form of Equation 
(H.14) and (H.15) yields 
















or 


1 

1 


a o 


"l 

5 

3.3480 


-°1- 


_2.6057 


The solution of this sst of squations givss ■ 1.4493 and a Q ■ 0.4493 so 
that tha polynomial can ba writtan as 

Z 2 - 1.4493 Z + 0.4493 - 0 . 

Tha roots of this polynomial ara Z^ ■ 1.0 and Z 2 • 0.4493, giving polas of 
8^^ • 0.0 and s 2 ■ 4.0003. Tha arror in ia dua to truncation arror. Nota 
that tha constralnad pola was ratumad axactly. 

H.3 M athod 2 

Anothar approach which can ba usad to constrain Prony's mathod to usa 
information about known polas is tha modification of tha polynomial root 
finding routina MULLER. Onca tha unconstrainad Prony's mathod calculatas 
tha coafficiants of tha polynomial (H.8) and if soma of tha roots of tha 
polynomial ara known a priori, than tha locations of thosa roots can ba 
passad to MULLER and it will not hava to saarch for thosa roots. Sines tha 
root finding routina la vary tlma consuming, tha knowlsdgs of tha location 
of any of tha roots will presumably save computation time. 

A possible flaw with this approach is that MULLER is forced to pre¬ 
sume roots of the polynomial when those exact roots may not ba contained in 
the polynomial. That is, if the polynomial has not been constrained to 
contain the known roots, as per Section H.2, then the known roots will 




not necessarily be contained in it. Forcing an unconstrained polynomial 
to have certain roots could grossly perturb the location of the other roots 
being sought. 

H.A Method 3 

The obvious solution to the flaw presented in Section H.3 is to use 
both the methods of H.2 and H.3 simultaneously. That 1st the polynomial 
is constrained to contain the known poles as outlined in Section H.2. Then 
the polynomial root finding routine MULLER can be modified to extract the 
known roots from the polynomial before it begins its search for the unknown 
roots. 

It is felt that this approach will give the best accuracy and will 
spaed up the calculations since the order of the polynomial is effectively 
reduced. 
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APPENDIX I 

REVERSING THE WAVEFORM IN TIME TO ELIMINATE 
EXTRANEOUS RESONANCES 


Reversing the waveform in time dapanda on cha atatiaticai properties 
of cha noiaa which do not changa whan Cha waveform ia ravaraad. Tha polaa 
chac raaulc from cha noiaa, cherafora, ara not alcarad by ravaraal in cima 
whlla cha crua polaa ara nagacad or flippad chrough cha origin. 

If cha noiaa laval ia high the noiaa polaa only approximately remain 
cha aame under Cima ravaraal and tha true poles only approximately reflect 
chrough the origin. If Cima ravaraal ia attempted for curve-fitting Prony's 
method it is found that all polaa flip precisely. Therefore, this method 
ia not affective for curve-fitting Prony's method. 

(Curve-fitting Prony'a method uaaa tha solution to cha Inhomogeneous system 
QX a q, in tha notation of Volume 1, where Q is a square matrix.) 

Tha waveform of Figure 1.1 waa uaad in a numerical example of the time-!* 
ravaraal method. Least-squares Prony's method waa applied to tha waveform. 
Estimates at the S-Plana polaa ware found using the inhomogeneous solution 
which is defined in Volume I, Section 2 as 

\ " [Q^r 1 <fq. 

The waveform consists of 100 samples and was corrupted with uncorrelated, 
Gauaslan-diatributed noise with a standard deviation of 0.1. Figure 1.2 
displays the poles obtained from the forward and reversed waveforms. The 
dimensions of the (NXn)-dimensional matrix Q in this example ara M*76 
and r„"24. Tha estimates of the true poles are quite accurate. Figure 1.3 
displays the poles obtained for tha same waveform but different matrix 
dimensions. The matrix dimensions are M*60 and n*40. In this case, the 
faxtraneoua poles do not all remain in the left-half plane which is 
attributed to using an overly square matrix. 
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If the matrix Is too "thin", of M/n »1, the estimates of the true 
poles become Inaccurate. If the matrix is too "fat", or M/n “ 1, the 
extraneous poles do not all remain in the left-half plane. The matrix 
shape where M/n « .1 appears to be the best shape for good estimates and 
keeping the extraneous poles in the left-half plana. 
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Figure X~3 Effect of an overly-square data matrix. 








Therefore, it appears that if we choose M/n * 3, then time 
reversal can be used to distinguish between the true poles and the 
extraneous poles in relatively high noise levels. 

The above two examples show that there is a trend for the 
noise poles in the least-squares Prony's method to occupy a higher 
frequency portion of the S-plane and to be highly damped. This trend 
was studied by letting the least-squares Prony's method operate on a 
waveform consisting of only Gaussian-dlstributed uncorrelated noise - 
no signal. The poles resulting from this example are highly damped and 
evenly distributed in the z-plane approximately around a circle within 
the unit circle. 


This behavior can be explained by the fact that the polynomial 

coefficients, except which is set to one, all tend to zero for the 

least-squares method operating on uncorrelated noise. The polynomial 

N 

then tends toward z ■ e where e tends to zero. The roots of this 
polynomial, z^, have magnitudes 




It then follows that the poles should be highly damped, since 
|z^| + 0 if |c| *0, and that the poles are evenly distributed about the 
z-plane. However, for curve-fitting Prony's method the a^, i ■ 0, ..., 
N-l do not ten 1 to zero and indeed this phenomenon is not observed in 
tests using cuive-fitting Prony’s method on all noise. 
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This behavior of the noise poles in least-squares Prony seems to 
remain approximately the same even if the waveform is not entirely noise. 

The trend that is seen is that the noise poles occupy the higher frequencies, 
are highly damped and evenly distributed between the higher frequencies; 
while the true poles are approximately at their uncorrupted locations 
and occupy the lower frequencies. That is, it appears as though the noise 
poles are "crowded" away from the lower frequencies or, perhaps more 
accurately, the lower frequency noise poles become lower frequency signal 
poles when the waveform is no longer entirely noise. 





